setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207")
filea=read.csv("all.FDR.sig.at.least.one.add.direction.same.diff.csv",header=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.csv("at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",header=T)
file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3 =file1[file1$BF_in_DC>10,]

enhancer1=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E071.bed",head=F,sep="\t")
enhancer2=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E074.bed",head=F,sep="\t")
enhancer3=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E068.bed",head=F,sep="\t")
enhancer4=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E069.bed",head=F,sep="\t")
enhancer5=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E072.bed",head=F,sep="\t")
enhancer6=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E067.bed",head=F,sep="\t")
enhancer7=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E073.bed",head=F,sep="\t")
enhancer8=read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/enhancer_p10/regions_enh_E070.bed",head=F,sep="\t")

enhancer=rbind(enhancer1,enhancer2,enhancer3,enhancer4,enhancer5,enhancer6,enhancer7,enhancer8)
enhancer =enhancer[!duplicated(enhancer$V4),]

promoter1 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E067.bed",head=F,sep="\t")
promoter2 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E068.bed",head=F,sep="\t")
promoter3 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E069.bed",head=F,sep="\t")
promoter4 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E070.bed",head=F,sep="\t")
promoter5 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E071.bed",head=F,sep="\t")
promoter6 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E072.bed",head=F,sep="\t")
promoter7 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E073.bed",head=F,sep="\t")
promoter8 =read.table("/mnt/md1200/7/zhaocunyou/promoter_enhancer/promoter_p10/regions_prom_E074.bed",head=F,sep="\t")
promoter =rbind(promoter1,promoter2,promoter3,promoter4,promoter5,promoter6,promoter7,promoter8)
promoter = promoter[!duplicated(promoter$V4),]

data = file2
data$snp.location=""
for(i in 1:dim(data)[1]){
chr=as.character(data[i,]$Chr)
pos=as.numeric(data[i,]$Start)
tmp1 =enhancer[enhancer$V1==chr,]
tmp1 =tmp1[tmp1$V2<=pos & tmp1$V3>=pos,]
tmp2 =promoter[promoter$V1==chr,]
tmp2 =tmp2[tmp2$V2<=pos & tmp2$V3>=pos,]
data[i,]$snp.location=ifelse(dim(tmp1)[1]>=1,paste0(data[i,]$snp.location,"enhancer"),paste0(data[i,]$snp.location,""))
data[i,]$snp.location=ifelse(dim(tmp2)[1]>=1,paste0(data[i,]$snp.location,"promoter"),paste0(data[i,]$snp.location,""))
}
#data1=data.frame(data$unitID,data$snp.location)
write.csv(data1,"../20210120/807.psy.ASH.add.enh.promtr.csv",quote=F,row.names=F)

#
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210120.H3k.analysis")

###进行一致性的分析
deepsea=read.table("./807.analysis/tmp/b67429bd-f7a1-49cc-af92-13c38d2125a8_807psyASH.vcf_FEATURE_zscore.tsv",head=T,sep="\t")
coln=names(deepsea)
braincell=c("Brain_Angular_Gyrus","Brain_Anterior_Caudate","Brain_Cingulate_Gyrus","Brain_Germinal_Matrix","Brain_Hippocampus_Middle","Brain_Inferior_Temporal_Lobe","Brain_Mid_Frontal_Lobe","Brain_Substantia_Nigra")

file1=read.csv("./807.analysis/consider.all.beta0.807.site.alt.up.down.statis.csv",header=T)
file2=read.csv("./807.analysis/807.psyASH.add.enh.promtr.csv",header=T)
file2=data.frame(unitID=file2$unitID,id=file2$id,snp.location=file2$snp.location)

file=merge(file1,file2,by="id")
locations=as.character(unique(file$snp.location))

for(lc in locations){
ASH=file[file$snp.location==lc,]#按位置进行分类
for(i in 1:length(braincell)){
braincn=coln[grep(coln,pattern = braincell[i])]
braincn=braincn[grep(braincn,pattern = "H3K")]
deepseq=deepsea[,c("chrom","pos",braincn)]
result=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
H3k=unlist(strsplit(braincn,"\\."))[seq(2,3*length(braincn),3)]
colnames(result)=H3k
frt=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
colnames(frt)=H3k
deepseq$id=paste0(deepseq$chrom,":",deepseq$pos)
deepseq=merge(deepseq,ASH,by="id")
for(j in 1:length(braincn)){
result[,j]=ifelse(deepseq[,3+j]<0,"down","up")
frt[,j]=ifelse(result[,j]==deepseq$group,"same","opposite")
}
rt_statis=data.frame(matrix(NA,length(braincn),4))
row.names(rt_statis)=H3k
colnames(rt_statis)=c("same","opposite","same_ratio","pvalue")
for(k in 1:length(braincn)){
rt_statis[k,1]=table(frt[,k]=="same")[2]
rt_statis[k,2]=table(frt[,k]=="opposite")[2]
rt_statis[k,3]=rt_statis[k,1]/(rt_statis[k,1]+rt_statis[k,2])
rt_statis[k,4]=binom.test(rt_statis[k,1],rt_statis[k,1]+rt_statis[k,2],p=0.5)$p.value
}
fn=paste0("./807.analysis/consider.all.beta0.",braincell[i],".",lc,"_statis.csv")
write.csv(rt_statis,fn,quote=F,row.names = T)
}
	#算均值然后求差异
sel1=paste0("./807.analysis/consider.all.beta0.",braincell,".",lc,"_statis.csv")
i=1
rt=read.csv(sel1[i],header = T)[,1:3]
for (i in 2:length(braincell)) {
  f1=read.csv(sel1[i],header = T)[,1:3]
  rt=rbind(rt,f1)
}

H3k.group=unique(rt$X)
col_names=c("H3k.group","same.mean","opposite.mean","same.ratio","P.value")
result=data.frame(matrix(NA,1,ncol=5))
names(result)=col_names
result=result[-1,]
for (i in 1:length(H3k.group)) {
  resultmp=data.frame(matrix(NA,1,ncol=5))
  names(resultmp)=col_names
  tmp=rt[rt$X==H3k.group[i],]
  resultmp[1,1]=tmp[1,1]
  resultmp[1,2]=round(mean(tmp[,2]),digits = 0)
  resultmp[1,3]=round(mean(tmp[,3]),digits = 0)
  resultmp[1,4]=resultmp[1,2]/(resultmp[1,2]+resultmp[1,3])
  resultmp[1,5]=binom.test(resultmp[1,2],(resultmp[1,2]+resultmp[1,3]),p=0.5)$p.value
  result=rbind(result,resultmp)
}
fn2=paste0("./807.analysis/807.ASH.",lc,".statis.csv")
write.csv(result,fn2,quote=F,row.names = F)
}

									#8544 ASH
deepsea=read.table("./8544.analysis/tmp/dd9a8167-9c49-4bf0-a3f7-15df0bbde2b2_8544.vcf_FEATURE_zscore.tsv",head=T,sep="\t")
coln=names(deepsea)
braincell=c("Brain_Angular_Gyrus","Brain_Anterior_Caudate","Brain_Cingulate_Gyrus","Brain_Germinal_Matrix","Brain_Hippocampus_Middle","Brain_Inferior_Temporal_Lobe","Brain_Mid_Frontal_Lobe","Brain_Substantia_Nigra")

file1=read.csv("./8544.analysis/consider.all.beta0.8544.site.alt.up.down.statis.csv",header=T)
file2=read.csv("./8544.analysis/8544.ASH.add.enh.promtr.csv",header=T)
file2=data.frame(unitID=file2$unitID,id=file2$id,snp.location=file2$snp.location)

file=merge(file1,file2,by="id")
locations=as.character(unique(file$snp.location))

for(lc in locations){
ASH=file[file$snp.location==lc,]#按位置进行分类
for(i in 1:length(braincell)){
braincn=coln[grep(coln,pattern = braincell[i])]
braincn=braincn[grep(braincn,pattern = "H3K")]
deepseq=deepsea[,c("chrom","pos",braincn)]
result=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
H3k=unlist(strsplit(braincn,"\\."))[seq(2,3*length(braincn),3)]
colnames(result)=H3k
frt=data.frame(matrix(NA,dim(ASH)[1],length(braincn)))
colnames(frt)=H3k
deepseq$id=paste0(deepseq$chrom,":",deepseq$pos)
deepseq=merge(deepseq,ASH,by="id")
for(j in 1:length(braincn)){
result[,j]=ifelse(deepseq[,3+j]<0,"down","up")
frt[,j]=ifelse(result[,j]==deepseq$group,"same","opposite")
}
rt_statis=data.frame(matrix(NA,length(braincn),4))
row.names(rt_statis)=H3k
colnames(rt_statis)=c("same","opposite","same_ratio","pvalue")
for(k in 1:length(braincn)){
rt_statis[k,1]=table(frt[,k]=="same")[2]
rt_statis[k,2]=table(frt[,k]=="opposite")[2]
rt_statis[k,3]=rt_statis[k,1]/(rt_statis[k,1]+rt_statis[k,2])
rt_statis[k,4]=binom.test(rt_statis[k,1],rt_statis[k,1]+rt_statis[k,2],p=0.5)$p.value
}
fn=paste0("./8544.analysis/consider.all.beta0.",braincell[i],".",lc,"_statis.csv")
write.csv(rt_statis,fn,quote=F,row.names = T)
}
	#算均值然后求差异
sel1=paste0("./8544.analysis/consider.all.beta0.",braincell,".",lc,"_statis.csv")
i=1
rt=read.csv(sel1[i],header = T)[,1:3]
for (i in 2:length(braincell)) {
  f1=read.csv(sel1[i],header = T)[,1:3]
  rt=rbind(rt,f1)
}

H3k.group=unique(rt$X)
col_names=c("H3k.group","same.mean","opposite.mean","same.ratio","P.value")
result=data.frame(matrix(NA,1,ncol=5))
names(result)=col_names
result=result[-1,]
for (i in 1:length(H3k.group)) {
  resultmp=data.frame(matrix(NA,1,ncol=5))
  names(resultmp)=col_names
  tmp=rt[rt$X==H3k.group[i],]
  resultmp[1,1]=tmp[1,1]
  resultmp[1,2]=round(mean(tmp[,2]),digits = 0)
  resultmp[1,3]=round(mean(tmp[,3]),digits = 0)
  resultmp[1,4]=resultmp[1,2]/(resultmp[1,2]+resultmp[1,3])
  resultmp[1,5]=binom.test(resultmp[1,2],(resultmp[1,2]+resultmp[1,3]),p=0.5)$p.value
  result=rbind(result,resultmp)
}
fn2=paste0("./8544.analysis/8544.ASH.",lc,".statis.csv")
write.csv(result,fn2,quote=F,row.names = F)
}